Identification of Immunological Biomarkers of Atopic Dermatitis by Integrated Analysis to Determine Molecular Targets for Diagnosis and Therapy

Introduction

Atopic dermatitis (AD) is a common chronic inflammatory skin disease, characterized by recurrent eczematous lesions and intense itching. AD is ranked as the leading cause of the global burden from skin inflammatory disease,¹ affecting approximately 10–30% of children and 2–10% of adults worldwide.² AD patients were divided into infantile AD, childhood AD and adult AD. Adult AD patients contain three forms, the persistent form of childhood-onset of AD, the relapsing form of childhood-onset of AD, and the adult-onset AD. Nummular eczema-like phenotype and prurigo nodularis-like pattern appeared to be more related to adult-onset AD, while lichenified/exudative flexural dermatitis was more common in childhood-onset AD.³ However, the mode of inheritance and genes involved in AD remain unclear. Currently, diagnosis relies exclusively on clinical features due to the lack of specific laboratory or histological findings.

In addition, conventional therapies were unsatisfactory for AD patients, due to the long-term side effects of glucocorticoid and immunosuppressant drugs. Recently, numerous monoclonal antibodies such as nemolizumab, blocking IL-31 receptor, dupilumab, blocking IL-4 receptor, lebrikizumab targeting IL-13 receptor have been widely used and proven effective in the treatment of AD. And new topical molecules such as JAK inhibitor, tofacitinib, phosphodiesterase 4 inhibitors, crisaborole and apremilast, and selective JAK1 and JAK2 inhibitor, baricitinib can also improve the clinical outcome of AD patients.⁴ Therefore, exploring specific targeting drugs opened a promising new era for AD treatment.

Immune cells, particularly dendritic cells (DCs) and CD4+ memory T cells, play key roles in AD progression. Because the immune response in AD involves functionally distinct cell types, it is necessary to assess and determine immune-infiltration changes specific to the condition to develop novel immunotherapeutic drugs that target various immune cells. As an analytical tool, CIBERSORT can be used to estimate the proportion of different immune cell types based on RNA-sequencing data.⁵ CIBERSORT contains 22 cell types, including B cells, T cells, monocytes, eosinophils, macrophages, natural killer cells, neutrophils, plasma cells, dendritic cells (DCs), and mast cells.⁶ The tool has been widely used to assess immune cell infiltrations in malignant tumors such as colorectal cancer and renal cell carcinoma.^7,8 In the AD skin tissue, DCs require cytokines, cytokine receptors, and T cell receptor (TCR) signals to recognize and present antigens to naïve T cells. These processes induce adaptive immunity to cause skin inflammation and barrier dysfunction. Furthermore, the impaired skin barrier facilitates infiltration of foreign antigens (eg, food allergens) as well as activation of pattern recognition and innate immune receptors. Consequently, IgE- and FCεRI-bearing dermal DCs are activated and migrate to the regional lymph nodes, inducing Th2 differentiation and B-cell IgE skewing. Specific T-cell subsets are recruited to the lesioned skin tissue, producing cytokines and chemokines. Both Th2 and Th22 (IL-22) cytokines inhibit epidermal differentiation and lipid synthesis, facilitating Staphylococcus aureus colonization and accelerating the so-called “AD march”⁹ Although characterized by Th2 immune responses, AD is now considered a heterogeneous disease, with additional activation of Th22, Th17/IL23, and Th1 cytokine pathways.¹⁰ However, the immune landscape of AD has not been entirely revealed. Therefore, we used CIBERSORT to explore effective immunological biomarkers for the diagnosis and treatment of AD.

In particular, we explored the immune landscapes of lesioned (LS) tissue and non-lesioned (NL) tissue from patients with AD and used CIBERSORT to determine their differences. Common differentially expressed genes (DEGs), hub proteins and their transcription factors (TFs), microRNAs (miRNAs), and long non-coding RNAs (lncRNAs) were analyzed using the GSE140684 and GSE99802 microarray datasets. Based on comprehensive bioinformatics analysis, we further aimed to investigate the molecular mechanisms underlying AD development and identify potential AD biomarkers (Supplementary Figure 1).

Materials and Methods

Microarray Data Processing and Identification of DEGs

From the GEO database, we obtained five expression profile datasets of AD, GSE140684, GSE99802, GSE58558, GSE120721and GSE36842. 50 samples of GSE140684,¹¹ including 29 from LS tissue and 21 from NL tissue, 112 samples GSE99802,¹² including 59 from LS tissue and 53 from NL tissue, 35 samples GSE58558,¹³ including 18 from LS tissue and 17 from NL tissue, 30 samples GSE120721,¹⁴ including 15 from LS tissue and 15 from NL tissue, 16 samples GSE36842,¹⁵ including 8 from LS tissue and 8 from NL tissue were selected in this study.

The original expression matrix was processed and analyzed by R software. DEGs were screened out by using the function lmFit and eBayes in the Limma package.¹⁶ Genes with adjusted P < 0.05 and log2FC > 1, were considered as DEGs of GSE140684, while genes with adjusted P < 0.05 and log2FC > 0.5 were considered as DEGs of GSE99802. TB tools was used to present Heatmap and perform Venn analysis to identify common DEGs from GSE99802 and GSE140684.¹⁷

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) is a widely used computational method for genome-wide expression matrix.¹⁸ In the present study, genetic information of GSE99802 and GSE140684 were uploaded to GSEA software (version 4.0.3). GSEA was performed to identify biological process (BP) Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways connected with AD. GSEA analysis was conducted according to default parameters. P < 0.05 was considered significant.

The enrichment analysis of DEGs was performed by using the clusterProfile package of R software and visualized by Hiplot (https://hiplot.com.cn/).¹⁹ KEGG pathways and GO terms with P<0.05 were selected.

Protein-Protein Interaction (PPI) Analysis

In order to provide insights into the mechanisms of AD, DEGs were uploaded to STRING (http://string-db.org) to construct PPI network.²⁰ Cytoscape (version 3.8.0), an open-source bioinformatics software platform was used to visualize molecular interaction networks and select key nodes.²¹ CytoHubba, a Cytoscape plugin, ranked the top 25 genes as hub genes, according to Maximal Clique Centrality (MCC).²²

TFs-mRNA Interaction Analysis

JASPAR database was used to identify regulatory TFs that influence the hub genes at a transcriptional level.²³ Hub gene-TF regulatory network was constructed and visualized by the Cytoscape software. The correlation between these TFs and hub genes was evaluated by Pearson correlation coefficient.

miRNA-mRNA and miRNAs-lncRNAs Interaction Analysis

The miRNA-mRNA regulatory network was constructed via miRTarBase v8.0,²⁴ including experimentally supported miRNA-gene interactions, and was further processed with Cytoscape. Then, miRNAs targeting more than 3 genes were selected. Based on StarBase 2.0,²⁵ the upstream molecules lncRNAs of the selected miRNAs were predicted with CLIP-Data ≥7, and were further analyzed by Cytoscape.

Cross-Validation of Candidate Biomolecules

DisGeNET is one of the largest publicly available platforms of genes and variants associated with human disease.²⁶ ImmPort is an open human immunology database for translational and clinical research.²⁷ In order to identify key genes and validate our workflow, we cross-checked the identified hub genes and TFs with DisGeNET and ImmPort. Then, we cross-checked the identified miRNAs by HMDD [22].

Drug-Hub Gene Interaction Analysis

To explore potential drugs for the treatment of AD, DrugBank database²⁸ was used to identify drug-hub gene interaction via NetworkAnalyst.²⁹ The identified target network was visualized by Cytoscape.

Prediction of Protein Subcellular Localization

We applied Cell-PLoc-2 to predict the subcellular localization of the proteins encoded by hub genes.³⁰ Cell-PLoc-2 is an improved online tool for predicting subcellular localization of proteins based on amino acid sequence.

Evaluation of Immune Cell Infiltration and Correlation Analysis of Immune Cell Infiltration Types and Hub Genes

To evaluate the abundance of infiltrating immune cells in LS skin and NL skin of AD, we used CIBERSORT gene signature file LM22, including neutrophils, dendritic cells (DCs), natural killer cells, macrophages, B cells, T cells, and subdivided resting and activated immune cell subtypes.⁶ CIBERSORT is based on linear support vector regression. Gene expression datasets GSE14064 and GSE99802 were processed and uploaded to CIBERSORT webservice portal (http://cibersort.stanford.edu/). P-value <0.05 was set as cut-off criteria. Differentially expressed immune cells between the LS group and NL group were obtained. Pearson correlation coefficient analysis was conducted and visualized by Hiplot in order to explore the relationship between these different subpopulations of immune cells and hub genes.

Construction of Least Absolute Shrinkage and Selection Operator (LASSO) Model and Receiver Operating Characteristic (ROC) Analysis

To select the best features for high-dimensional data, we applied LASSO with strong predictive value and low correlation. The expression profile of GES99802 hub genes were extracted to construct LASSO model by glmnet (https://CRAN.R-project.org/package=glmnet), so as to distinguish LS skin from NL skin of AD. A model index was established using the regression coefficients from the LASSO analysis to weight the expression value of the hub genes. Then, Hiplot was used to conduct and perform ROC curve analysis in GSE140684, GSE99802, GSE58558, GSE120721, and GSE36842.³¹

Results

Identification of DEGs

Based on the sample information and data matrix, 551 and 445 genes were differentially expressed in LS and NL samples of GSE140684 (adj.P< 0.05, logFC>1) and GSE99802 (adj.P<0.05, logFC>0.5), respectively. The 20 most up- and down-regulated gene heatmaps were visualized in volcano plots (Figure 1A–D). 238 DEGs common to both GSE140684 and GSE99802 were found (Figure 1E).

Figure 1 Identification of DEGs in microarray datasets GSE140684 and GSE99802. (A) Heat map of top 20 down-regulated and up-regulated DEGs of GSE140684 (B) Heat map of top 20 down-regulated and up-regulated DEGs of GSE99802. (C) Volcano plot of the distributions of all DEGs of LS and NL samples of GSE140684 (adj. P < 0.05, logFC >1). (D) Volcano plot of the distributions of all DEGs of LS and NL samples of GSE99802 (adj. P < 0.05, logFC >0.5). Blue dots represent significantly down-regulated genes and red dots represent significantly up-regulated genes in LS samples. (E) Venn plot of common DEGs identified between two data sets. Abbreviations: FC, fold-change; DEGs, differentially expressed genes; NL, non-lesioned; LS, lesioned.

DEGs from LS Samples Were Mainly Enriched in the Cytokine-Cytokine Receptor Interaction Pathway and Immune Response

According to GSEA-based KEGG analysis, gene sets associated with the cytokine-cytokine receptor interaction and chemokine signaling significantly enriched in LS samples of the GSE140684 and GSE99802 datasets (Figure 2). In addition, GSEA-based GO biological analysis revealed that the most strongly enriched terms were regulation of immune responses, DCs, and CD4+ T cells (Table 1).

Table 1 Significantly Enriched GSEA-Based GO Terms of Biological Process Related to Immune Cells and Responses with P-value <0.05 Were Screened Out

Figure 2 GSEA-based KEGG analysis. (A) GSEA-based KEGG analysis of GSE140684. (B) GSEA-based KEGG analysis of GSE99802. Abbreviations: GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

KEGG analysis of the 238 DEGs revealed numerous genes enriched in immune-related pathways (cytokine-cytokine receptors interaction, viral protein interaction with cytokine and cytokine receptor, cell adhesion molecules, chemokine signaling, IL-17 signaling, T cell receptor signaling, and NF-kappa signaling) (Figure 3A). Among these, the cytokine-cytokine receptor pathway was the most significantly enriched pathway with 28 associated genes, including TNFSF8, CCR5, CCR7, CSF2RA, FLT3, CXCL1, CXCL2, IL2RG, IL4R, IL7R, CXCL8, IL13RA2, LTB, PRLR, CCL8, CCL17, CCL19, CCL20, CCL22, CXCL1, IL18RAP, IL27RA, CXCR6, IL24, IL37, IL21R, IL26, and IL36G. BP analysis revealed that the DEGs were mainly related to T cell activation, positive regulation of cell-cell adhesion, leukocyte cell-cell adhesion, positive regulation of leukocyte cell-cell adhesion, neutrophil migration, positive regulation of T cell activation, and neutrophil chemotaxis (Figure 3B).

Figure 3 Enrichment analysis of the common DEGs identified from GSE140684 and GSE99802. (A) Top seven enriched KEGG terms. (B) Top seven enriched GO biological process terms. Abbreviations: DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Proteomic Signatures in AD

To identify the core genes among the DEGs identified between LS and NL samples, we constructed the PPI network and obtained 198 nodes plus 762 edges (Figure 4A). A high MCC indicates that a given protein is clustered in the center of the network and plays an essential role. We identified 25 hub genes in the PPI, including SELL and SLAMF1; we then used the MCC to select genes for further analysis (Figure 4B). The hub genes were mainly related to the immune system, including cytokines, cytokine receptors, TCR signaling molecules, selectins, and integrins (Table 2).

Table 2 Summary of Hub Proteins Identified from Protein-Protein Interactions Analysis of Encoded Differentially Expressed Genes in Atopic Dermatitis

Figure 4 Identification of hub genes. (A) PPI network of the DEGs in Atopic dermatitis. The nodes indicate the DEGs and the edges indicate the interactions between two proteins. Medium confidence score was used for the construction of PPI networks. (B) Top 25 hub genes with MCC by CytoHubba plugin. Abbreviations: PPI, protein-protein interaction; MCC, maximum correlation criterion.

Competing Endogenous RNA (ceRNA) Signatures

We constructed the miRNA-mRNA regulatory network (Figure 5A) and selected seven miRNAs (gene counts ≥3) to predict interacting lncRNAs (Figure 5B). Among these miRNAs, we further investigated the roles of hsa-mir-124-3p, hsa-mir-21-5p, and hsa-mir-155-5p in AD (Table 3).

Table 3 Summary of miRNAs (Genes Count ≥3) of Regulatory Biomolecules of the Hub Genes in Atopic Dermatitis Identified from Hub Genes-miRNA s Interactions

Figure 5 Construction of ceRNA regulatory networks. (A) miRNA-mRNA regulatory network. Medium confidence score was used for the construction of regulatory networks. Genes are colored in blue, and node size is adjusted according to number of targeted miRNAs; miRNAs are colored in purple; miRNAs targeting more than two genes simultaneously are colored in red. (B) miRNA-LncRNA regulatory network. MiRNA with genes count ≥ 3 were screened to predict LncRNA with CLIP-Data ≥7. Abbreviations: ceRNA, competing endogenous RNA; miRNA, micro-RNA; LncRNA, long non-coding RNA.

Correlation Between TFs and Hub Genes

The TF-hub gene regulatory networks (Figure 6A) identified potential roles of STAT3, JUN, PPARG, and RELA in AD (Table 4). The TF that targeted the highest number of hub genes was FOXC1, which showed significantly downregulated expression in both datasets. We then analyzed correlations between FOXC1and its targeted hub genes (CCL20, CXCL1, IL2RG, CCR7, CTLA4, CD2, SLAMF1, LCK, LCP2, and ICOS) in the GSE90882 and GSE140684 datasets. The FOXC1 expression level was negatively correlated with the production of chemokines (CCL20, CXCL1), cytokine receptors (IL2RG, CCR7), and TCR signaling molecules (CTLA4, CD2, SLAMF1, LCK, LCP2, ICOS) (Figure 6B). These significant correlations in two independent datasets further indicated that these hub genes and FOXC1 play an indispensable role in AD development.

Table 4 Summary of TFs (Genes Count ≥5) of Regulatory Biomolecules of the Hub Genes in Atopic Dermatitis Identified from Hub Genes-TFs Interactions

Figure 6 Construction of TF regulatory networks and correlation analysis. (A) TF-mRNA regulatory networks. Medium confidence score was used for the construction of regulatory networks. (B) Correlation analysis between TF FOXC1 and targeted genes. Abbreviation: TF, transcription factor.

Drug Analysis

Seven hub genes (CCR5, CXCL8, MMP9, LCK, ITGAL, ITK, IL2RG) had interactions with 42 drugs for possible treatment (Figure 7A).

Figure 7 Construction of protein-drug interaction network and subcellular localization. (A) Protein-drug interaction network of hub genes. (B) The distribution and percentages of the subcellular localization of the proteins encoded by hub genes.

Employing Cell-PLoc-2, we identified the subcellular localization of proteins encoded by the 25 hub genes of AD. The proportions of hub proteins in different subcellular locations are displayed in Figure 7B. Notably, the hub nodes LCP2 and ITK displayed cytoplasm localization; CD2, ITGAX, and IKZF1 were localized in the nucleus; SELL, CTLA4, CD28, SLAMF1, CD3E, CD5, IL2RG, CCR5, and CD274 were localized in plasma membrane; and IL7R, ITGAL, CXCL8, CCL19, ICOS, CCL20, and CXCL1 were localized in the extracellular space.

Correlation Analysis Between Hub Gene Expression and Immune Cell Infiltration

Based on the CIBERSORT output results, we analyzed 42 (84%) samples from GSE99802, while excluding GSE140684 entirely because it contained 21 (42%) non-significant (p > 0.05) samples, which prevented performing an accurate analysis. The different populations of infiltrating immune cells in the LS and NL groups are presented using a violin plot in Figure 8A. Three types of immune cells (resting CD4+ memory T cells, activated CD4+ memory T cells, and DCs) were significantly differentially expressed in both datasets, and were specifically upregulated in the LS tissue.

Figure 8 Violin plot of immune cell composition between LS and NL tissue and correlation analysis between immune cells and hub genes. (A) The violin plot indicates the composition of 22 immune cells between LS and NL tissue in GSE140684 with CIBERSORT p < 0.05 for all eligible samples. The blue violin plot indicates NL tissue, and the red violin plot represents LS tissue. (B) Correlation analysis between differentially expressed immune cells and hub genes. Abbreviations: NL, non-lesioned; LS, lesioned.

To further explore the relationship between hub gene expression and immune cell infiltration, we analyzed whether the expressions of hub genes were correlated with the levels of resting and activated CD4+ memory T cells as well as activated DCs. The levels of activated and resting CD4+ memory T cells were moderately correlated with the expressions of chemokine-encoding genes (eg, CXCL8, CXCL1, CCL20) and cytokine receptor-encoding genes (eg, IL7R, IL2R, CCR5). CCL19 and CCR7 expression levels were positively correlated with the levels of activated and resting CD4+ memory T cells, respectively. Additionally, expression levels of TCR signaling molecules (CTLA4, CD3E, SLAMF1, LCK, LCP2, ICOS, and ITK, but not CD274) were significantly correlated with the levels of activated and resting CD4+ memory T cells (Figure 8B). We also found a positive correlation between the levels of activated DCs and the expressions of all hub genes, except for CCL20.

Exploring Candidate Biomarkers by ROC Curves and LASSO Regression

The LASSO regression model was firstly constructed for the GSE99802 dataset (with the largest number of samples) to search for an optimal linear combination of hub genes that can be used to diagnose AD (Figure 9A). Six genes were selected with non-zero regression coefficients and lambda.min = 0.0383513. The LASSO regression model was established as follows: index = SELL*0.343074243+IL7R*0.004220025+CXCL1*0.165991590+CCR7*0.416155475+CCR5*0.007626045+CCL19*0.050059214. Each hub gene and the LASSO model was evaluated according to the ROC curve in five cohorts (GSE99802, GSE140684, GSE58558, GSE120721, and GSE36842) (Figure 9B). The AUC of the LASSO model was in the range of 0.75–1.00. CD274, GZMB, and CXCL1 showed particularly high performance, with AUC values above 0.9 in GSE36842, GSE120721, and GSE120721, respectively, indicating their potential as biomarkers of AD.

Figure 9 Construction of LASSO regression model and ROC curves of hub genes in five cohorts. (A) The left plot indicates binomial deviance of different numbers of variables revealed by the LASSO regression model for GSE99802. The red dots represent the value of binomial deviance; the grey lines represent the SE; the vertical dotted lines represent optimal values by the minimum criteria and 1-SE criteria. “Lambda” is the tuning parameter. (B) The ROC curves of LASSO regression model (SELL *0.343074243+IL7R*0.004220025+CXCL1*0.165991590+CCR7* 0.416155475 +CCR5*0.007626045+ CCL19*0.050059214) and top 5 genes in 5 cohorts (GSE99802, GSE140684, GSE58558, GSE120721, GSE36842). Abbreviations: AUC, area under the curve; ROC, receiver operating characteristic; SE, Standard error.

Discussion

The multifactorial background of AD occasionally causes therapeutic failure, necessitating individualized treatment and new solutions. One recommendation is to find a method that can distinguish patients with AD according to various characteristics such as disease severity, biomarkers, and immunological polarization (ie, phenotypes, immunotypes).¹⁰ This strategy would allow for specific treatments to different forms of AD. Toward this end, in this study, we employed CIBERSORT to identify novel biomarkers and immunotypes that can be used to distinguish patients with AD. Enrichment analysis revealed AD-associated genes were mainly related to the immune response and cytokine interactions. Moreover, our results suggest that inhibition of the TF FOXC1 might upregulate hub genes, including those encoding cytokines, cytokine receptors, and TCR signaling molecules, to promote antigen presentation and activation of DCs. Surface receptors and TCR signaling molecules of CD4+ T cells can recognize antigens presented by DCs and initiate their differentiation into Th2 cells, causing local inflammation, epidermal dysfunction, and tissue remodeling and fibrosis, thus contributing to AD pathogenesis. In addition, the LASSO regression model of hub genes can be used to identify diagnostic biomarkers.

Our exploration of AD-related biological processes and pathways began with GSEA-based GO and KEGG pathway analyses. We found that AD-affected tissues were enriched in immune cells (CD4+ T cells and DCs), cytokine-cytokine receptors, as well as chemokine and TCR signaling pathways. It is well recognized that AD is a Th2/Th22-dominant disease, regulated by DC-induced T-cell polarization and recruitment of specific T-cell subsets by chemokines.³² Th17 and Th1 skewing also plays a role in causing AD.³³ Th2 cells release cytokines such as IL-4, IL-5, and IL-13, leading to elevated IgE production that increases skin inflammation and aggravates the skin barrier defect in AD.³⁴ These results were consistent with previous research and emphasized the important role of Th17 and Th1 polarization, apart from Th2/Th22 axis.

MicroRNAs are small endogenous non-coding RNAs that bind to the 3′ untranslated region of genes and regulate gene expression through degrading or inhibiting the translation of target genes.³⁵ Our results identified several hub genes, miRNAs, confirmed by previous research. AD is characterized by impaired skin barrier and chronic inflammation.³⁶ MMP-9, as broad-spectrum protease, can break down a range of matrix and other proteins, causing substantial tissue damage and barrier dysfunction of AD pathology.³⁷ Elevated level of CCR5 in AD induced chronic inflammation, impairing the skin of AD.³⁸ Additionally, by targeting the TF PU.1, microRNA-155 regulates DC-associated Th2 responses and controls allergic inflammation in mice.³⁹ Our study indicated these molecules were worthy of further exploration and might act as promising therapeutic targets.

Then, we used CIBERSORT to estimate the composition of immune cells and further correlate these cells with hub genes. Among the hub genes, cytokine-, cytokine-receptor‒, and TCR-encoding genes were positively associated with infiltration of CD4+ memory T cells and activated DCs. Activated DCs can migrate to the regional lymph nodes, where they induce Th2 production and subsequent B-cell IgE switching. We found that activated DCs were positively correlated with the hub genes CCR5, CXCL1, and CXCL8. It was reported that CCL19 binds to CCR7 of DCs in the AD LS skin, leading to DC activation and Th2 polarization.⁴⁰ Additionally, DCs can promote cancer progression via CXCL1 and aggravate CXCL8 in chronic obstructive pulmonary disease.^41,42 Together with these findings, our results might suggest that activated DCs could release CXCL1 and CXCL8 to recruit neutrophils and cause skin inflammation in AD. Besides, activated DCs are positively correlated with TCR signaling molecules, suggesting that DCs upregulate these molecules and promote T cell differentiation, thereby impairing epidermal integrity and causing cutaneous inflammation.

Following antigen stimulation, T cells persist as effector memory T cells and central memory T cells in the lymph nodes and peripheral tissues.⁴³ Antigen recognition by TCR and interaction of inflammatory cytokines such as IL-2 with cytokine receptors are required for the generation and maintenance of tissue resident memory (TRM) T cells. IL7R and IL2RG can bind to IL-2, promoting the generation of CD4+ memory T cells.⁴⁴ Additionally, CCR5+ TRM cells can differentiate into Th17 and Th1 cells.⁴⁵ CCL19 promotes the differentiation of CCR7+ CD4+ memory T cells into Th2 subsets,⁴⁶ whereas CCL20 induces CD4+ memory T cell differentiation into Th17 cells, leading to increased CXCL1 and CXCL8 production.⁴⁷ In addition, CD4+ memory T cell numbers are significantly related to TCR molecules (except for CD274); antigen recognition of TCRs is essential for CD4+ T cell differentiation into memory subsets.⁴⁸ Our findings firstly suggest that CD4+ memory T cells might cause epidermal barrier dysfunction and contribute to AD by impairing intercellular cohesion and integrity, and disturbing keratinocyte differentiation.

Correlation analysis between differentially expressed TF FOXC1 and immune cell-related hub genes revealed important TFs that regulate hub gene expression, thus enhancing our understanding of AD pathogenesis. FOXC1, a member of the forkhead box TF family, is necessary for human keratinocyte terminal differentiation.⁴⁹ Previous research showed that CCL20 might induce the differentiation of CD4+ memory T cells into Th17 cells, attracting neutrophils, DCs, and more T lymphocytes to perpetuate skin inflammation. CCL19 might activate DCs and CD4+ memory T cells, leading to Th2 polarization in AD lesions. Our results suggest that FOXC1 might play a role in AD by upregulating the production of cytokines, cytokine receptors, and TCR signals. However, these conclusions require further investigation and confirmation.

We applied ROC analysis to evaluate single hub genes or LASSO linear models. The model consists of SELL, IL7R, CXCL1, CCR7, CCR5, and CCL19; however, three single genes (GZMB, CD274, and CXCL1) showed great diagnostic value in skin tissue.

Our study possesses several strengths, but also some limitations. We applied several novel analytic bioinformatic approaches on micro-array data of different datasets. The results were cross validated by DisGeNet and HMDD. A LASSO model was constructed based on one dataset and validated by four external datasets. However, we did not carry out experiments to validate the results. In order to better illustrate the key link between LS and NLl skin of AD patients, we did not include healthy controls in this study and ignore the basal genetic difference between individual. Therefore, a larger cohort, including healthy controls and more functional experiments are required to confirm the data.

Conclusion

In summary, our study employed next generation sequencing method to explore the differentially expressed genes of AD. Cibersort and LASSO model were applied to evaluate immune landscapes and diagnostic values of these genes in AD, respectively. A series of novel potential biomarkers, including immunological hub genes, TFs, miRNAs, and lncRNAs, were identified as candidate biomarkers and promising therapeutic targets in the prediction and diagnosis of AD. And the specific molecular function of these molecules needs to be further investigated.